!--------------------------------------------!
!        A. Gazizov, LNGS, 27.12.2010        !
!  included in project ZN on 27.11.12, LNGS  !
!          coefficients on IR/UV of          !
!  $ \gamma + p \rightarrow e^+ e^- p $ and  !
!  $ \gamma + N \rightarrow N' + X $         !
!          divided by H_0                    !
!--------------------------------------------!

module ircoeff
!   use iruv, only: nEIR, EIR, WIR
!   use approx1x
   implicit none
   
   integer, parameter:: neeir = 94   &
                      , nNXir = 99   &
					  , ncir  = 2    &  ! p->p or p->n; ...
					  , nmir  = 6    &  ! [1,2..6]=(P,B,D,Bs,Ds,Dss)
                      , nzir  = 51

   real(8):: ln10, lnE0                 ! ln(10) and norm. energy E0=1E18 eV
   real(8), dimension(neeir)        :: lnEeeir
   real(8), dimension(nNXir)        :: lnENXir
   real(8), dimension(nzir )        :: zzir
   real(8), dimension(2,neeir,nzir) :: beeir 
   real(8), dimension(ncir,nmir,nNXir,nzir):: PIR  ! pp/pn; 1,2..6; .. .. 

   integer:: ir_ready=0

!------------------
     contains
!------------------

!---------------------------------------------
Subroutine ReadIRCoeff
   implicit none
   integer:: i, j, k, ic, im

   ln10 = dlog(1d1)
   lnE0 = 1.8d1*ln10	! E0 in eV

! array of z
   do j = 1,nzir
      zzir(j) = 1d-1*(j-1)
   enddo

! Energies of protons
   do i = 1,nNXir  
      if(i.le.53) then 
	     lnENXir(i) = (15.8d0+0.05d0*(i-1))*ln10
      else
         lnENXir(i) = (18.5d0+0.1d0*(i-54))*ln10
	  endif
   enddo

! Coefficients of pair-production on EBL
   open(1,file='DataIR/beepIR.dat')
      read(1,'(103E15.8)') (lnEeeir(i),((beeIR(k,i,j),k=1,2)  &
                           , j=1,nzir),i=1,neeir)
   close(1)
   print*, ' File DataIR/beepIR.dat is read'

! Coefficients of nuclear photoproduction on EBL
   open(1,file='DataIR/PIRtkn.dat')
      read(1,'(99E15.8)') ((((PIR(ic,im,i,j),i=1,nNXir)  &
	                        , j=1,nzir),im=1,nmir),ic=1,ncir)
   close(1)
   print*, ' File DataIR/PIRtkn.dat is read'

   ir_ready = 1

end Subroutine ReadIRCoeff
!------------------------------


!---------------------------------------------------
! Get comoving $\beta_{ee}(E,z)/H_0$  
! and $H_0^{-1} d\beta_{ee}(E,z)/dlnE$ 
! on EBL of T. Kneiske et al. best fit.
! To be multiplied by $(1+z)^3/Redshift(1+z)$
! r = ln(E/E0), where E_0 = 10^{18} eV 

Subroutine GetbeeIR(r,z,bir)
   use msimsl, only: DQD2VL
   implicit none
   LOGICAL:: CHECK
   real(8),intent( in):: r, z
   real(8),dimension(neeir,nzir):: tmp
   real(8),intent(out),dimension(2):: bir

   real(8):: rr, ztmp
   integer:: k

   if(ir_ready.eq.0) call ReadIRCoeff 

   if(z.gt.5d0) then
      ztmp=5d0
   else
      ztmp = z
   endif

   rr = r + lnE0    ! ln(E) in eV

   if(rr.lt.lnEeeir(1)) then
      bir  = 0d0
	  return
   endif

   if(rr.gt.lnEeeir(neeir)) rr = lnEeeir(neeir)
   do k = 1,2      
      tmp = beeir(k,:,:)
      bir(k) = DQD2VL(rr,ztmp,neeir           &
             , lnEeeir,nzir,zzir,tmp,neeir,CHECK)
   enddo
   bir(2) = bir(1)*bir(2)

end Subroutine GetbeeIR
!--------------------------------------------------

!---------------------------------------------------
! Get comoving 
! P_pp, B_pp, D_pp, B'_pp, D'_pp, D''_pp [1.. 6]
! P_pn, B_pn ,D_pn, B'_pn, D'_pn, D''_pn [7..12]
! all $P_{pi}(E,z)/H_0$  
! and $B'_{pi}=H_0^{-1} dB_{pi}(E,z)/dlnE$ 
! on EBL of T. Kneiske et al. best fit.
! To be multiplied by $(1+z)^3/Redshift(1+z)$
! r = ln(E/E0), where E_0 = 10^{18} eV 

Subroutine GetPIR(r,z,P)
   use msimsl, only: DQD2VL
   implicit none
   LOGICAL:: CHECK
   real(8),intent( in):: r, z
   real(8),intent(out),dimension(ncir,nmir):: P
   real(8),dimension(nNXir,nzir):: tmp

   real(8):: rr, ztmp
   integer:: ic, im

   if(ir_ready.eq.0) call ReadIRCoeff 

   if(z.gt.5d0) then
      ztmp=5d0
   else
      ztmp = z
   endif

   rr = r + lnE0    ! ln(E) in eV

   if(rr.lt.lnEeeir(1)) then
      P = 0d0
      return
   endif

   if(rr.gt.lnEeeir(neeir)) rr = lnEeeir(neeir)

   do ic = 1,ncir
      do im = 1,nmir
         tmp = PIR(ic,im,:,:)
         P(ic,im) = DQD2VL(rr,ztmp,nNXir           &
                  , lnENXir,nzir,zzir,tmp,nNXir,CHECK)
      enddo ! im
   enddo ! ic

end Subroutine GetPIR
!--------------------------------------------------


end module ircoeff